Load required packages

library(tidyverse)
library(phyloseq)
library(speedyseq)
library(ggrepel)
library(ampvis2)
library(plotly)
library(microbiome)
options(getClass.msg=FALSE) # https://github.com/epurdom/clusterExperiment/issues/66
#this fixes an error message that pops up because the class 'Annotated' is defined in two different packages

Load functions from Github

'%!in%' <- function(x,y)!('%in%'(x,y))

source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_taxa_tests.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_normalisation.R")
## Loading required package: scales
## 
## Attaching package: 'scales'
## The following object is masked from 'package:microbiome':
## 
##     alpha
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
## Loading required package: reshape2
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_alpha.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_beta.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_heatmap.R")

Load physeq object

ps = "data/processed/physeq_update_11_1_21.RDS"

ps %>% 
  here::here() %>%
  readRDS() %>%
  phyloseq_get_strains_fast() %>%
  phyloseq_remove_chloro_mitho() -> physeq
## Joining, by = "ASV"
physeq %>% 
  subset_samples(Experiment == "Continuous") %>% 
  subset_samples(Paul %!in% c("Paul")) %>%
  subset_samples(Reactor != "IR2") -> ps_PolyFermS

We will be analysing only the PolyFermS samples here so take a subset of the physeq object.

physeq %>% 
  subset_samples(Experiment == "Continuous") %>% 
  subset_samples(Paul %!in% c("Paul")) %>%
  subset_samples(Reactor != "IR2") -> ps_polyFermS

sample_data(ps_polyFermS)$Reactor <- fct_relevel(sample_data(ps_polyFermS)$Reactor, "IR1", "CR", "TR1", "TR2","TR3", "TR4", "TR5", "TR6") 

sample_data(ps_polyFermS)$Treatment <- fct_relevel(sample_data(ps_polyFermS)$Treatment, "UNTREATED",  "CTX+HV292.1", "CTX","HV292.1","VAN+CCUG59168", "VAN",  "CCUG59168") 

sample_data(ps_polyFermS)$Reactor_Treatment <- fct_relevel(sample_data(ps_polyFermS)$Reactor_Treatment, "IR1_UNTREATED","CR_UNTREATED", "CR_CTX", "CR_VAN", "TR1_CTX+HV292.1","TR2_CTX", "TR3_HV292.1", "TR5_VAN+CCUG59168", "TR4_VAN", "TR6_CCUG59168") 

ps_polyFermS %>% 
  rarefy_even_depth(sample.size = 4576,
                    rngseed = 123) -> ps_rare
## `set.seed(123)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(123); .Random.seed` for the full vector
## ...
## 16 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## CR-40-S166ETR1-30-S178ETR1-42-S194ETR2-30-S195IR1-40-S197
## ...
## 50OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
ps_rare %>%
  subset_samples(Enrichment == "NotEnriched") %>%
  subset_samples(Day_from_Inoculum >= 22 & Day_from_Inoculum <= 38) -> ps_polyFermS_rare

ps_polyFermS_rare
## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 296 taxa and 72 samples ]:
## sample_data() Sample Data:        [ 72 samples by 50 sample variables ]:
## tax_table()   Taxonomy Table:     [ 296 taxa by 8 taxonomic ranks ]:
## phy_tree()    Phylogenetic Tree:  [ 296 tips and 294 internal nodes ]:
## refseq()      DNAStringSet:       [ 296 reference sequences ]
## taxa are rows

Alpha div:

Compute alpha-div metrics:

ps_polyFermS_rare %>%
  phyloseq_alphas(phylo = TRUE) -> alpha_df
measures = c("Observed", "diversity_shannon", "evenness_pielou", "PD","MNTD", "SES.MPD" ,"bNTI")


alpha_df %>%
  plot_alphas(measure = measures,
             x_group = "Reactor_Treatment",
             colour_group = "Reactor",
             fill_group = "Reactor",
              shape_group = NULL,
              facet_group = "Reactor_Treatment",
            test_group = "Reactor_Treatment",
            test_group_2 = NULL) -> out

out$plot + 
  facet_null() + 
  facet_grid(alphadiversiy ~ ., scales = "free") + 
  ggpubr::rotate_x_text(60) +
  scale_fill_viridis_d() +
  scale_color_viridis_d()

out$stat %>%
  dplyr::filter(signif == "SIGN") %>%
  DT::datatable()

Beta-div:

Compute beta-div metrics:

ps_polyFermS_rare  %>%
  phyloseq_plot_bdiv(bdiv_list,
                     m = "CoDa",
                     seed = 123) -> coda
  
coda$PCA$layers[[1]] = NULL

coda$PCA + geom_point(size=2,
                   aes(color = Reactor_Treatment, 
                       fill = Reactor_Treatment,
                       shape = NULL,
                       alpha = Day_from_Inoculum)) + 
  theme_light() +
  # geom_path(arrow = arrow(type = "open", angle = 30, length = unit(0.15, "inches")),
              # size = 0.05, linetype = "dashed", inherit.aes = TRUE, aes(group=Reactor_Treatment, color = Reactor_Treatment), show.legend = FALSE) +
  scale_alpha_continuous(range=c( 0.9, 0.3)) + scale_color_viridis_d(na.value = "black") +
  scale_fill_viridis_d(na.value = "black") + 
  # scale_shape_manual(values = c(8, 21, 22, 23, 24, 16, 15, 18, 17)) + 
  theme_classic() +
    labs(col=NULL, fill = NULL, shape = NULL) + guides(shape=FALSE) -> p1

p1

p1 %>%
  ggplotly() -> p1ply

p1ply
ps_polyFermS_rare %>%
  phyloseq_compute_bdiv(norm = "pc",
                        phylo = TRUE,
                        seed = 123) -> bdiv_list
## Loading required package: ape
## Loading required package: GUniFrac
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
## 
## Attaching package: 'vegan'
## The following object is masked from 'package:microbiome':
## 
##     diversity
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
## 
##     count
bdiv_list$coda <- coda$physeq_clr %>% phyloseq::distance(method = "euclidean")
phyloseq_plot_bdiv(dlist = bdiv_list, # list of distance computed from a phyloseq object
                   ps_rare = ps_polyFermS_rare, # phyloseq object
                   m = "PCoA", # PCoA or NMDS
                   seed = 123, # for reproducibility
                   axis1 = 1, # axis to plot
                   axis2 = 2) -> plot_list
## [1] "bray"
## [1] "bjaccard"
## [1] "wjaccard"
## [1] "uunifrac"
## [1] "wunifrac"
## [1] "d_0"
## [1] "d_0.5"
## [1] "coda"
plot_list$bray = NULL; plot_list$d_0 = NULL; plot_list$d_0.5 = NULL

plot_list %>%
  phyloseq_plot_ordinations_facet(color_group = "Reactor_Treatment",
                                 shape_group = NULL) +
  scale_color_viridis_d(na.value = "black") +
  scale_fill_viridis_d(na.value = "black") 

plot_list %>%
  phyloseq_ordinations_expl_var()
## New names:
## * .id -> .id...1
## * V1 -> V1...2
## * .id -> .id...3
## * V1 -> V1...4
##    .id...1           V1...2  .id...3           V1...4
## 1 bjaccard   Axis.1   [17%] bjaccard Axis.2   [14.7%]
## 2 wjaccard Axis.1   [22.8%] wjaccard Axis.2   [15.8%]
## 3 uunifrac Axis.1   [20.2%] uunifrac Axis.2   [15.6%]
## 4 wunifrac Axis.1   [67.8%] wunifrac Axis.2   [14.4%]
## 5     coda Axis.1   [19.2%]     coda Axis.2   [11.9%]
ps_polyFermS_rare %>%
  phyloseq_distance_boxplot(bdiv_list$coda ,
                            d = "Reactor_Treatment") -> dist_box


dist_box$plot

  lapply(
  bdiv_list,
  FUN = physeq_pairwise_permanovas,
  physeq = ps_polyFermS_rare,
  compare_header = "Reactor_Treatment",
  n_perm = 999,
  strat = FALSE
) %>%
  bind_rows(.id = "Distance") %>%
  mutate_if(is.numeric, round, 3) %>%
  # filter(! terms %in% (c("Residuals", "Total"))) %>%
  DT::datatable()

Taxa Heatmap:

ps_polyFermS_rare %>%
  # subset_samples(Experiment == "CCUG59168" |  Experiment ==  "HV292.1" | Experiment ==  "Cecum" |  Experiment ==  "Continuous" & Reactor %in% c("IR1", "CR")) %>%
  # rarefy_even_depth(sample.size = 2574,rngseed = 123) %>%
  phyloseq_ampvis_heatmap(transform = "percent",
                          group_by = "Day_from_Inoculum",
                          facet_by = c("Enrichment", "Phase", "Reactor", "Treatment", "Experiment", "Reactor_Treatment" ),
                          tax_aggregate = "Class",
                          tax_add = NULL,
                          ntax  = 60) -> p
## Warning: There are only 7 taxa, showing all
p + facet_grid( ~ Reactor_Treatment + Phase , scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2

ps_polyFermS_rare %>%
  # subset_samples(Experiment == "CCUG59168" |  Experiment ==  "HV292.1" | Experiment ==  "Cecum" |  Experiment ==  "Continuous" & Reactor %in% c("IR1", "CR")) %>%
  # rarefy_even_depth(sample.size = 2574,rngseed = 123) %>%
  phyloseq_ampvis_heatmap(transform = "percent",
                          group_by = "Day_from_Inoculum",
                          facet_by = c("Enrichment", "Phase", "Reactor", "Treatment", "Experiment", "Reactor_Treatment" ),
                          tax_aggregate = "Genus",
                          tax_add = NULL,
                          ntax  = 60) -> p

p + facet_grid( ~ Reactor_Treatment + Phase , scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2

ps_polyFermS_rare %>%
  # subset_samples(Experiment == "CCUG59168" |  Experiment ==  "HV292.1" | Experiment ==  "Cecum" |  Experiment ==  "Continuous" & Reactor %in% c("IR1", "CR")) %>%
  # rarefy_even_depth(sample.size = 2574,rngseed = 123) %>%
  phyloseq_ampvis_heatmap(transform = "percent",
                          group_by = "Day_from_Inoculum",
                          facet_by = c("Enrichment", "Phase", "Reactor", "Treatment", "Experiment", "Reactor_Treatment" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 60) -> p

p + facet_grid( ~ Reactor_Treatment + Phase , scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2

paste0(here::here(),
       "/data/processed/",
       "stab",
       "_",
       format(Sys.time(), "%Y%b%d")
       ,".RData") %>% save.image()
# load("/Users/fconstan/Projects/EZe/ASV/260420.RData")
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] GUniFrac_1.1         matrixStats_0.57.0   vegan_2.5-7         
##  [4] lattice_0.20-41      permute_0.9-5        ape_5.4-1           
##  [7] metagMisc_0.0.4      reshape2_1.4.4       scales_1.1.1        
## [10] microbiome_1.10.0    plotly_4.9.3         ampvis2_2.6.4       
## [13] ggrepel_0.9.0        speedyseq_0.5.3.9001 phyloseq_1.32.0     
## [16] forcats_0.5.0        stringr_1.4.0        dplyr_1.0.3         
## [19] purrr_0.3.4          readr_1.4.0          tidyr_1.1.2         
## [22] tibble_3.0.5         ggplot2_3.3.3        tidyverse_1.3.0.9000
## 
## loaded via a namespace (and not attached):
##   [1] Rtsne_0.15          colorspace_2.0-0    ggsignif_0.6.0     
##   [4] rio_0.5.16          ellipsis_0.3.1      rprojroot_2.0.2    
##   [7] XVector_0.28.0      fs_1.5.0            rstudioapi_0.13    
##  [10] farver_2.0.3        ggpubr_0.4.0        DT_0.17            
##  [13] ggnet_0.1.0         fansi_0.4.2         lubridate_1.7.9.2  
##  [16] xml2_1.3.2          codetools_0.2-18    splines_4.0.3      
##  [19] doParallel_1.0.16   knitr_1.30          ade4_1.7-16        
##  [22] jsonlite_1.7.2      broom_0.7.3         cluster_2.1.0      
##  [25] dbplyr_2.0.0        compiler_4.0.3      httr_1.4.2         
##  [28] backports_1.2.1     assertthat_0.2.1    Matrix_1.3-2       
##  [31] lazyeval_0.2.2      cli_2.2.0           htmltools_0.5.1    
##  [34] prettyunits_1.1.1   tools_4.0.3         igraph_1.2.6       
##  [37] gtable_0.3.0        glue_1.4.2          Rcpp_1.0.6         
##  [40] carData_3.0-4       Biobase_2.48.0      cellranger_1.1.0   
##  [43] zCompositions_1.3.4 vctrs_0.3.6         Biostrings_2.56.0  
##  [46] multtest_2.44.0     nlme_3.1-151        crosstalk_1.1.1    
##  [49] iterators_1.0.13    xfun_0.20           network_1.16.1     
##  [52] openxlsx_4.2.3      rvest_0.3.6         lifecycle_0.2.0    
##  [55] rstatix_0.6.0       zlibbioc_1.34.0     MASS_7.3-53        
##  [58] hms_0.5.3           parallel_4.0.3      biomformat_1.16.0  
##  [61] rhdf5_2.32.2        RColorBrewer_1.1-2  curl_4.3           
##  [64] yaml_2.2.1          NADA_1.6-1.1        stringi_1.5.3      
##  [67] S4Vectors_0.26.1    foreach_1.5.1       BiocGenerics_0.34.0
##  [70] zip_2.1.1           truncnorm_1.0-8     rlang_0.4.10       
##  [73] pkgconfig_2.0.3     evaluate_0.14       Rhdf5lib_1.10.0    
##  [76] labeling_0.4.2      patchwork_1.1.1     htmlwidgets_1.5.3  
##  [79] tidyselect_1.1.0    here_1.0.1          plyr_1.8.6         
##  [82] magrittr_2.0.1      R6_2.5.0            IRanges_2.22.2     
##  [85] generics_0.1.0      picante_1.8.2       DBI_1.1.0          
##  [88] foreign_0.8-81      pillar_1.4.7        haven_2.3.1        
##  [91] withr_2.4.0         mgcv_1.8-33         abind_1.4-5        
##  [94] survival_3.2-7      car_3.0-10          modelr_0.1.8       
##  [97] crayon_1.3.4        PhyloMeasures_2.1   rmarkdown_2.6      
## [100] progress_1.2.2      grid_4.0.3          readxl_1.3.1       
## [103] data.table_1.13.6   reprex_0.3.0        digest_0.6.27      
## [106] stats4_4.0.3        munsell_0.5.0       viridisLite_0.3.0